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Fourier Techniques for Very Long Astrophysical Time Series Analysis 

Scott M. Ransom^'^, Stephen S. Eikenberry'^, and John Middleditch* 

ABSTRACT 

o , 

. In this paper, we present an assortment of both standard and advanced Fourier techniques 

that are useful in the analysis of astrophysical time series of very long duration — where the 
pi ' observation time is much greater than the time resolution of the individual data points. We 

<^ . begin by reviewing the operational characteristics of Fourier transforms (FTs) of time series data, 

including power spectral statistics, discussing some of the differences between analyses of binned 
data, sampled data, and event data, and briefly discuss algorithms for calculating discrete Fourier 
transforms (DFTs) of very long time series. We then discuss the response of DFTs to periodic 
signals, and present techniques to recover Fourier amplitude "lost" during simple traditional 
analyses if the periodicities change frequency during the observation. These techniques include 
, Fourier interpolation which allows us to correct the response for signals that occur between Fourier 

' frequency bins. We then present techniques for estimating additional signal properties such as 

the signal's centroid and duration in time, the first and second derivatives of the frequency, the 
P\J \ pulsed fraction, and an overall estimate of the significance of a detection. Finally, we present a 

' recipe for a basic but thorough Fourier analysis of a time series for well-behaved pulsations. 

Subject headings: pulsars — methods: data analysis 

o 

, 1. Introduction 

The analysis of time series data is an important tool in many areas of astrophysics, including research 
' involving white dwarfs, black holes, and neutron stars. In the study of neutron stars, time series analysis has 

■ had particular importance for pulsar research due to the high coherence of pulsar periodicities. While many 

techniques can be used to investigate the properties of these periodic signals, the computational efficiency 
of the Fast Fourier Transform (FFT) makes Fourier analysis the preferred approach for many applications 
(see e.g. Burns & Clark 1969). The literature contains literally hundreds of descriptions of various aspects 
of Fourier analysis, most of which deal with signal detection using the power spectrum. 

Because of this concentration on power spectra, many researchers discard a wealth of information 
provided by the Fourier phases. Techniques which use this phase information exist that can provide insight 
into many useful signal properties. While many of these techniques have been known for some time (see e.g. 
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Middleditch 1976), few have appeared in textbooks or refereed journals, and fewer still have been presented 
with any sort of derivation or insight into their assumptions and/or limitations. 

A second, more practical problem with most astronomical Fourier analysis is its concentration on short 
time series. We define "short" to mean either that the full time series of binned or sampled data can fit into 
the core memory of your computer (N < 10^ points), or for data consisting of events (such as photon arrival 
times in X-ray astronomy) , that the time resolution (dt) of each event makes up a non-negligible fraction of the 
total time duration (T) of the data {T /dt < 10^). Billion point FFTs have been used successfully in the past 
(e.g. Anderson 1992; Mattox et al. 1996; Middleditch et al. 2000), but each required the use of state-of-the-art 
supercomputing facilities. Today, such analyses are possible using clusters of workstations or even individual 
desktop machines. Many projects, such as pulsar searches of globular clusters, astero- or hclioseismological 
observations, and gravitational wave experiments, require extremely large Fourier transforms in order to 
make the highest sensitivity (i.e. coherent) searches for pulsations, and to extract the maximum amount of 
information from signals found in these searches. 

It is our goal in this paper to present some useful Fourier analysis techniques, that for various reasons are 
used only rarely when working with long time series. Most of our examples pertain to pulsar searches of very 
long time series, but the methods can be used in the Fourier analysis of virtually any coherent periodicity. 
This paper will briefly discuss the properties of the DFT, its response to periodic signals and noise, and 
methods for its computation. We will discuss methods for interpolating Fourier amplitudes, estimating 
a signal's duration and centroid in time, accurately determining its frequency and frequency derivative, 
correcting for changing pulsation frequency during an observation, and estimating the phase- or amplitude 
modulation of a signal. Such techniques allow the detection of signals whose frequencies or amplitudes change 
significantly during an observation — for instance, due to orbital motion or pulsar spin-down. 

2. Discrete Fourier Transform 

In order to present more advanced Fourier techniques later, we first review some fundamentals of the 
DFT and its most common implementation, the FFT. Since thorough discussions of the Fourier transform 
in both its continuous and discrete variants exist in the literature (e.g. Bracewell 1999), we will mention only 
a few topics of particular relevance to astrophysical data analysis, closely following Middleditch (1976). 

2.1. Introduction to the DFT 

The A;*'' element of the Discrete Fourier Transform (DFT), of a uniformly-spaced time series nj (j = 
0, 1, . . . , A'' — 1) is defined as 

N-l 

Ak=Yl e-^-'^'^/^, (1) 

where i = V— 1 and k is the Fourier frequency or wavenumber {k = 0,1, . . . , N — 1). For a time spacing dt 
between successive data elements, the frequency of the k^^ element is fk = k/(N dt) = k/T, where T is the 
total time duration of the sequence being transformed. This frequency spacing of 1/T, for evenly sampled. 
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un-padded and non-windowed data, is often called the Independent Fourier Spacing (IFS)^ , and defines the 
finest frequency resolution available while maintaining completely independent Fourier amplitudes. Fourier 
frequency N/2T is known as the Nyquist frequency. 

If we view the DFT summation in the complex plane, we see that it is a simple vector addition with 
each element rotated by —2'jTk/N from the previous element. If the nj have a constant value, the sum will 
form k regular iV/fc-sided polygons with each polygon returning near to the origin, and with the last one 
returning exactly to the origin. Therefore, the DFT of a constant data string will be identically zero for all 
fc > 0, and equal to the sum of the data elements for fc = (the "DC" frequency clement). 

For most astrophysical observations, the data points are real-valued. This property adds an important 
feature to Fourier transforms of these time series — they are symmetric about the Nyquist frequency, such 
that AjY-k = Al, where represents the complex conjugate of Ak- This symmetry allows us to calculate 
the ftdl DFT of a time series by computing amplitudes at only half of the normal Fourier frequencies, thereby 
speeding up computation of the DFT by nearly a factor of two (e.g.. Press ct al. 1992). 

When deriving properties and techniques based on DFTs, it is often both computationally and intuitively 
easier to work with a time-normalized data series, where T = 1 and u represents the fraction of the observation 
complete at any given instant (such that < m < 1). In this case, instead of working with frequencies /, or 
integral wavenumbcrs fc, we define r, which represents any real wavenumber (including non- integers). If the 
number of samples, A^, from our data source, n{u), is large, we can compute a continuous FT 



which is almost identical to our DFT when r = k, and produces very high accuracy estimates of the Fourier 
amplitudes at any frequency such that <C r <C N/2. We will use this approximation in many of our 
derivations. 



The Fast Fourier Transform is a family of well-known computer algorithms which quickly calculate a 
DFT in only 0(A^log2(iV)) operations as opposed to the 0{N^) operations of a brute-force DFT. FFTs, 
their computation, and their origins have been described in numerous articles and books over the last few 

decades (see Braccwell 1999, for an introduction). Therefore, wc will describe only a few special versions 
which have become generally known only recently and which are useful in the analysis of extremely long 
time series. 

High energy pulsar searches using photon counting systems (infrared, optical, x-ray, and gamma-ray 
detectors) or pointed radio telescope searches (for example, searching globular clusters) often utilize very 
high sampling rates (i.e. 10 — 50 kHz) and very long integration times of hours, days, or even weeks. These 
observations result in time series with hundreds of millions or even billions of data points. The subsequent 
FFTs are impossible to perform using conventional FFT algorithms unless the full time series fits into the 
core memory of the computer being used. Utilizing special "out-of-core" FFT methods, we can transform 



^The IFS is important when trying to determine the overall significance of a candidate from a search. The number of IFSs 
searched corresponds to the number of independent trials searched and should therefore be included in calculations that try to 
determine if a candidate was produced by noise. See Vaughan et al. (1994). 




(2) 



2.2. Computation of Very Long FFTs 
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such huge arrays using distributed memory parallel systems or with a single workstation and manageable 
numbers of passes through the data stored on external media. Most of these methods are based on the 
"four-step" FFT and/or special external media array permutation methods (Fraser 1976; Bailey 1990). 



2.2.1. The "Four-Step" FFT 

Closely following the derivation described in Sweet & Wilson (1995), we can think of our one-dimensional 
time series as a "C-likc" or row-ordcrcd, two-dimensional matrix of size N = NrNc, where Nc is the number 
of columns (i.e. the length of a row) and Nr the number of rows (i.e. the length of a column). Using this 
data decomposition, the FFT is computed by 1. FFTing the columns of the matrix, 2. multiplying the data 
by complex "twiddle factors", 3. FFTing the rows of the matrix, and 4- matrix transposing the result. If we 
define indices x = 0, 1, . . . , A^c — 1, y = 0, 1, . . . , Af^r — 1, I = 0,1, ... ,Nr — 1, and m = 0,1, ... ,Nc — 1, we 
can write our signal and its DFT amplitudes as 

n{x,y) =nj, j = Nc y + X (3) 

and 

A{l,m) ^ Ak, k ^ Nr m + l. (4) 
Substituting into the definition of the DFT (eqn. 1) and simplifying we get 



Ail,m)= 



x=0 



^-2nixi/N J2 n(a;,t/)e-2^'^'/^'- 



^-2-K\xm/Nc 



Notice that the bracketed terms are really Nc DFTs of length A^^ — FFTs of the matrix columns multiplied 
by the "twiddle factor" g-^'^'^'/^. We denote these column FFTs as 

Ac{x,l) = e-^'^'^'/^ "(^-y) e-'"'^'/^'-. (6) 

3/=0 

We now see that the outer summation is iV^ DFTs of length Nc — FFTs of the matrix rows composed of 
the Ac{x, I) terms. We denote the result of these transforms as 

Afc-l 

Ar{m,l) = Y ^c(x,0 e-2--™/JV^ (7) 

x=0 

To recover the full FFT in its normal order, we simply transpose Ar{m, I) 

A{l,rn) = A^r{rn,l). (8) 

The "four-step" algorithm only needs small portions of the data in memory at any one time. Unfor- 
tunately, the short FFTs are strided in memory^ by instead of being stored contiguously. This results in 
a significant inefficiency in today's cache-based processors and requires inter-node communications when 
distributing the short FFTs over many processors on parallel computer systems. These shortcomings can be 
overcome with the "six-step" algorithm (Bailey 1990). 



"Strided" means that the data is stored in a sequence of non-contiguous memory locations spaced by a constant amount of 
memory known as the "stride" . 
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2.2.2. The "Six-Step" FFT 

If the initial data set is transposed from a, NcX Nr matrix into a x Nc matrix, the strided column 
FFTs of length Nr become contiguous row FFTs. This memory locality facilitates the use of processor 
cache systems, greatly increasing memory response times, and allows independent calculation of the row 
FFTs by different processors in parallel systems. A second transpose operation after the row FFTs have 
been calculated counteracts the effects of the first transpose operation and makes the next set of row FFTs 
contiguous in memory as well. The FFT is finished with a final transpose operation which brings the data 
into normal order. 

These transpose operations are relatively efficient on distributed memory parallel systems with fast inter- 
node communications, and add little time to the overall FFT. In fact, for applications with very large N, the 
time required to move the data from external media to computer memory and vice versa tends to dominate 
the FFT time. The "six-step" algorithm has become the "standard" FFT algorithm for distributed memory 
systems where the serial nature and large number of short FFTs exploit parallel computation strengths. 

Unfortunately, if the data docs not all fit into the core memory of a workstation or parallel machine, 
the transposes become extremely slow operations since data must be read to and written from much slower 
external media. Fraser (1976) devised optimized methods for dealing with such data permutations on external 
media including the "two-pass" FFT algorithm. 

2.2.3. The "Two-Pass" Out-of-Core FFT 

Fraser (1976) and Bailey (1990) both describe how a very large data set may be Fourier transformed with 
only two read-write passes through externally stored data if a scratch area on the external media the same 
size as the input data set exists. The method uses the "four-step" algorithm with out-of-core transposes. 
These transpose algorithms allow blocked external media access and perform most of the transposition work 
in core memory. While external media access speeds and transfer rates are orders of magnitude slower than 
core-memory systems, the "two-pass" algorithm allows huge arrays to be transformed in manageable times. 

3. Fourier Transforms of Real Data 

3.1. Fourier Response to Noise 

If our data nj are composed of some constant value cj plus random real- valued noise fluctuations dj 
with zero mean, the transform terms become 

{cj + dj ) e-'-^^'^/^ = Cj e-^^'^^''' + dj 6-2-'^'=/^ . (9) 

The linear nature of the Fourier transform allows us to treat the DFT of the dj independently from the 
constant length steps, Cj. Since the complex phase factor for a given j and k is fixed, the direction of each 
element in the sum is nearly fixed. However, since the sign of the dj may be either positive or negative, the 
vector direction of the j**^ element may be reversed. Thus, the DFT of the dj creates a kind of random walk 
in the complex plane. 

The statistical properties of this random walk for DFTs of pure noise have been well studied (see, e.g. 
Blackman & Tukey 1959), and result in power spectra distributed according to an exponential distribution 
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(a distribution with 2 degrees of freedom) with average and standard deviation equal to N (d^^. If we 
normahze the powers by dividing by N the probabihty for a power P = in a single bin to equal 

or exceed a power P' by chance is'' 

Prob{P > P') = e-^' . (10) 

Similarly, if we sum m properly normalized powers, the probability for the summed power P^ to exceed a 
power P' is given by 

Prob{Pm>P')=Y.^—fe-'", (11) 

j=0 ^■ 

which is the probability for a distribution of 2m degrees of freedom to exceed 2P' . Such an incoherent 
(since no phase information is used) summation of powers is often useful when searching for signals suspected 
of having power in many harmonics (see §3.2.2). 

Proper normalization of the powers is essential for an accurate estimate of a signal's statistical signifi- 
cance or lack thereof. We often cannot normalize our power spectrum by simply dividing by N (dj), since 
frequency-dependent noise may be present throughout our power spectriim perhaps as a result of instru- 
mental, atmospheric, or astrophysical processes. Typically, these processes produce noise which increases in 
strength towards the low-frequency part of the spectrum and is correspondingly called red noise. 

Techniques to flatten or remove this "coloured" noise component from the power spectrum are described 
by Israel & Stella (1996), and usually involve dividing short portions of the power spectrum by the locally 
determined average power level, Piocai, such that 

P _ \Au\' ^, l^fcl' _ Pk 

N{d])- Plocal' Plocal- ^ ' 

As long as the number of averaged powers is small enough such that the power spectrum is roughly constant 
over the range in question, a white noise like power spectrum is produced with average and standard deviation 
of approximately one and an exponential distribution (eqn. 10). 

Since strong narrow-band signals near some frequency of interest will skew a local power average upwards 
(and correspondingly decrease the calculated significance of a signal detection), it is important to exclude 
such powers from the calculation of Piocai ■ A simple and effective way to accomplish this is by normalizing 
with a corrected local median power level instead of the local power average. An exponential distribution 
with unity mean and standard deviation has a median of ln(2). Therefore, dividing a section of raw powers 
by l/ln(2) times the local median value is theoretically equivalent to normalizing with the local mean, but 
has the advantage of being insensitive to high-power outliers in the spectrum. 

More advanced algorithms for the removal of "coloured" noise and power normalization do exist. A 

simple example involves fitting polynomial models to portions of the power spectrum and then dividing 
by them out. These methods work well for Fourier frequencies near zero where the assumption of roughly 
equivalent power levels for the local powers may be unwarranted. 

For the special case where the noise in our data is purely Poissonian (i.e. for binned photons in an 
optical or x-ray observation), we have = (rij). In this case, our properly normalized power for the A;*^ 



''This is difTerent than the probability for an actual signal to produce a power P > P' in the presence of noise (see §3.3). 
This difference is important in setting upper limits on the amplitudes of periodic signals as discussed in Vaughan et al. (1994). 



- 7- 



DFT element is 

P, - ^ - ^ (13) 

where nph = N (rij) is the sum of the nj (or the total number of photons for a photon-counting system), which 
is also equal to the "DC" frequency value of the FT (see §2.1). However, the same processes that caused 
the "coloured" noise discussed above can significantly alter this situation and require a power normalization 
based on local powers (see eqn. 12). 



3.2. Fourier Response to Periodic Signals 

One of the more useful properties of the FT for astronomical purposes is its response to periodic signals. 
Since all real periodic signals can be expanded into a series of sinusoids it is important to understand the 
FT response to a simple sine wave. 



3.2.1. Sinusoidal Signals 

If we now let our nj represent a sampled cosinusoid of amplitude a, phase ^, and frequency fr = r/T 
(where wavenumber r is an integer and fr an "integral frequency"), we can write 

nj = a cos{2-Kfrjdt + (j)) (14a) 
= a cos{2n jr /N + (/)) (14b) 

= ^ ^g27rijr/JV+i<^ _|_ g-27ry>/Ar-i<^^ _ ^j^^^-j 

From this expression, we see that the /c*''' element of the DFT is given by 

iV-l 

^ 3=0 
N-1 

= - ^2mj{k-r)/N+i(j> _^^-2-Kij{k+r)/N-i<p^ ^^r^^ 

and represents the summation of two vectors in the complex plane. For k ^ r, the first term traces out 
\k — r\ complete "rotations" (pseudo-polygons which start and end at the origin) in the complex plane (since 
fc — r is an integer), giving a net contribution of zero to the fc*^ DFT element. The second term traces out 
k + r complete rotations and once again contributes nothing to the fc*^ DFT element (since k and r are both 
positive and therefore k + r ^0. 

When k = r, however, the consecutive terms in the summation add coherently (i.e. in phase and therefore 
without rotation), since the rotation caused by the DFT exponential exactly cancels that from the signal 
exponential (the Fourier transform "derotates" the signal). As a result, each element of the sum is a step in 
the complex plane of magnitude a/2 in a direction parallel to the one set by the arbitrary initial phase of 
the signal, (f). For a cosimisoidal signal with integral frequency fr, the DFT will be uniformly zero, except 
in the r^^ frequency bin, where the response is 

Ar = ^ e'^. (16) 
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The Fourier response is more complicated for sinusoids witli non-integral frequencies (i.e. wavenumber r 
is a non- negative real number). The fc"^ DFT element is still given by eqn. 15b, but not all of the signal ends 
up in a single DFT bin A^. When k = [r] (where [r] is the nearest integer to r), the first term in eqn. 15b 
traces out a fraction (fc — r) of a complete rotation in the complex plane, while the second term traces out 
k + [r] complete rotations, plus an additional fractional rotation. 

When TV is large, these complete and fractional "rotations" can bo treated as circles and arcs respectively. 
Therefore, the first term of eqn. 15b results in a semi-circular arc of length Na/2 along the arc, while the 
second term produces a semi-circular arc of length a/2 (Nmodk + [r]) along the arc. The DFT response 
is simply a vector drawn from the origin to the end of the arc (see Fig. 1). Since virtually all astrophysical 
applications involve r » 1, where the first term dominates the response, we will ignore the second term in 
the rest of our analysis. 

The final response is a chord subtending 27r(A: — r) radians of a circle of radius Na/A'K{k — r). The 
equation for the length of a chord is 

^^ = 2|:sinf|V (17) 



e V2, 

where s is the arc length and 8 is the angle subtended by the chord. The curvature of the arc away from 
the signal's starting phase </> results in a phase change of e~™'^^~^\ Therefore, the DFT response and power 
are 



A, = e'^e-'-(^-'-)2 ,y . sm 



27r(fc - r) 



(18a) 



27r(fc - r) 
= ^^10 sinKfc- r)] 

2 iT{k-r) ^ ' 

= Ao e-^'^('=-'') sine [7r(fc - r)] (18c) 
Pk = l^fel^ = Po sinc'[7r(fc-r)], (18d) 

where Ao = Na/2 e"^ is the DFT response for an integral frequency signal (eqn. 16), Po is the corresponding 
Fourier power, and the sine function is defined as sinc(a;) = sin(a;)/a;. This result is easily confirmed by a 
direct integration of eqn. 2 where n{u) is equal to eqn. 14c with j/N — > u. 

The sine factor in eqn. 18c produces a loss of sensitivity for the standard FFT to most real-world signals 

(where r is not an integer). This effect, often called "scalloping" (e.g. Middleditch et al. 1993), is shown in 
Fig. 2, and causes a worst-case (when |fc — r| = 1/2) amplitude reduction of \Ak\ = 2 |Ao| /ti" — a nearly 
60% loss of signal power. On average, scalloping results in a ~ 23% loss of signal power (van der Klis 1989). 
It is important to remember, though, that this loss in sensitivity is due to the finite frequency resolution of 
the FFT algorithm rather than an intrinsic feature of the data itself. In §4.1 we discuss various methods to 
reduce or even eliminate this loss of sensitivity. 



3.2.2. Non- Sinusoidal Signals 

Many real- world periodic signals are not sinusoidal. Fortunately, we can expand all real- valued pulsations 
as a series of m sinusoidal components 

m 

Hj = ah cos{2TTj hr/N + (ph) (19a) 
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h=l 



where h = 1,2, ... ,m specifies the harmonic number (with h = 1 known as the "fundamental"), and Uh and 
(f)h represent the amplitude and phase of each component respectively. Due to the linear nature of the FT, 
we can treat the harmonics as independent sinusoidal signals. Each of these sinusoids produces a Fourier 
response equivalent to eqn. 18c, except that Ag becomes Ah = Nah/2 e'"^''. 

For nearly sinusoidal pulsations only the first few terms of eqn. 19a contain significant amplitudes, an- 

This results in a similarly small number of significant peaks in the corresponding power spectrum of the 
data. Low duty-cycle pulsations (i.e. those with a pulse that is short compared to the pulse period), such 
as most radio pulsars, on the other hand, have dozens of significant terms in their expansions and therefore 
harmonics in their power spectra. 

A useful pulsation model, particularly for radio and x-ray pulsars, can be constructed based on a modified 
von Mises distribution (MVMD) 

^KCOs(27rfrt-\-(f)) p — 

where <t <T is the instantaneous time, Iq is the modified Bessel function of zeroth order, and the shape 
parameter, k, determines the width of the function (e.g. Mardia & Zemroch 1975). In the limit as k — > 0, 
the MVMD becomes a sinusoid, while as k ^ oo, it becomes a Gaussian (see Fig. 3). The integral of the 
MVMD over a single pulse period is simply a, all of which is pulsed (i.e. the pulsed fraction is one). The 
full-width at half-maximum (FWHM) as a fraction of a pulse, is 

FWHMmvmd = TT"^ arccos {In [cosh(K)]} , (21) 

and the maximum value is 

2acosh(«;) , 

maxMVMD = :rT^ 7- (22) 

Io(k;) - e-*^ 

The FT of the MVMD can be computed in a particularly convenient form for harmonic analysis. 
According to Abramowitz & Stegun (1972, eqn. 9.6.34), we can expand the exponential in the MVMD 

as 

oo 

g«cos(x) = + 2 ^ lh{K) cos{hx), (23) 

h=l 

where I/i is the modified Bessel function of order h. When combined with the rest of the MVMD definition 

we have 

lo(K) -e " 

This expression is simply a "DC" term (since the integral over a pulse equals a) plus a series of independent 
cosinusoidal harmonics. After Fourier transforming (i.e. substituting into eqn. 2 with f^t ru), we are left 
with a series of harmonics of amplitude ANlh{K)/ [lo('v) — e"**] and phase hcj) at Fourier frequency hr. It is 
important to note that as k ^ oo and the pulse bec;omes narrower, the Fourier amplitudes of the low order 
harmonics are twice that of a sinusoid with the same pulsed fraction (see the dashed line in Fig. 4). This 
fact, along with the large number of harmonics that low duty-cycle pulsations generate, can significantly 
increase search sensitivities to such pulsations (see Fig. 5). 

Fig. 4 shows the approximate number of significant harmonics (meaning that a harmonic's amplitude 
is greater than one-half the amplitude of the fundamental) generated by an MVMD pulsation, as well as 
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a histogram of the duty-cycles of over 600 radio pulsars (the majority of which are from Taylor et al. 
1995). Most radio pulsars have duty-cycles < 5% corresponding to > 10 significant harmonics — assuming 
a sufficient data sample rate. 

3.3. Periodic Signals With Noise 

When a periodic signal is present in a noisy time series, a sum of m powers Pm, containing some amount 
of signal power Ps, is no longer described by a distribution with 2m degrees of freedom (sec §3.1). Groth 
(1975) calculated the expectation value and variance of P„ as {Pm) =m + P, and (^P^ - {P,nf) = m + 2P, 
respectfully. He also derived the exact probability density function for Pm which can be integrated to give 
the probability that Pm is greater than or equal to some power P', 

Proh{{Pm;Ps)>P')=e-i^'^'')Y: E (25) 

fc=0 J=0 ■'' 

When Ps ~Q this equation reduces to Eqn. 11. 

The fact that the probability density function for a signal plus noise is different from a distribution 
with 2m degrees of freedom is very important when trying to determine the sensitivity of a search for pulsa- 
tions or an upper limit to the amplitude of a periodic signal present in a time series. Vaughan et al. (1994) 
describe a procedure* for correctly determining search sensitivities and upper limits using the equations of 
Groth (1975). 

3.4. Photon Counting Data 

Since many of today's astronomical time series come from photon counting experiments, it is important 
to raise some of the issues particular to Fourier analysis of such data. If we can assume purely Poissonian 
statistics, a power spectrum of pure noise is flat, and can be normalized simply by dividing by the total 
number of photons in the data (the zeroth or "DC" frequency bin from the DFT — see §3.1). In addition 
to this difference in power spectrum normalization, the other points worth noting come from the fact that 
photon counting data is based on the measurement of events rather than the instantaneous sampling of a 
continuous process. 

One important issue that is beyond the scope of this paper is dead-time correction. Dead-time effects 
modify a detector's sensitivity to photons for some time after the detection of an earlier photon. These 
effects can cause complicated non- linear and frequency-dependent systematics during Fourier analysis. We 
refer the reader to Zhang et al. (1995) and references therein, for a thorough discussion of this topic. 

3.4-1- Binned vs. Sampled Data 

Many high-energy telescopes and detectors produce time series of binned photons rather than the sam- 
pled data produced by radio telescopes. Since binning essentially averages a periodic signal's instantaneous 



*Note that Vaughan et al. (1994) use a power normalization that is a factor of two higher than ours. 
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rate over the binning time (dt), it modifies the Fourier response to the signal. Binning removes phase informa- 
tion from the data and causes the Fourier response to sinusoidal pulsations to become frequency-dependent 
— resulting in decreased sensitivity at high frequencies (e.g. Middleditch 1976; Leahy et al. 1983). 

The frequency-dependent loss in Fourier amplitude due to binning is smc{nfrdt) or sinc(7rr/7V). The 
binned data Fourier response to a sinusoid is therefore eqn. 18c times this factor. This decrease in sensitivity 
corresponds to a loss in signal power of about 19.8% at half the Nyquist frequency and 59.5% at the 
Nyquist frequency itself. 

For Poissonian noise (i.e. from a photon-counting experiment that docs not introduce count-rate de- 
pendent systematics) , which is independent of the sampling interval, the Fourier response is flat over all 
frequencies. This is in contrast to a sinusoidal signal passing through the same system which suffers the 
frequency-dependent attenuation described above (Middleditch 1976). Such behavior is important when try- 
ing to estimate limits or amplitudes for pulsations in a time series (Vaughan et al. 1994). 



3.4-2. Low Count Rate Data 

The Fourier analysis of gamma-ray or x-ray observations often places us in a very unique regime — 
very long integration times (> 10'* s) with very low numbers of counts (< 10^ photons). In addition, due to 
visibility constraints based on the orbits of the telescopes, large fractions of the time between the first and 
last photons may be devoid of counts. 

Fourier analysis of such data can overwhelm present computational resources. For example, a 10^ s 
observation (about 11.6 days) with photon time-of-arrival (TOA) resolution of 10~^ s would require a 10 Gpt 
FFT for a full-resolution analysis. Such FFTs, while possible, are extremely difficult to compute unless very 
special and dedicated hardware resources are available. If this data contains only a small number of photons, 
however, we can exactly compute the DFT over any frequency range and to any frequency resolution using 
a brute-force implementation of the FT. 

If we treat each TOA as a sample of amplitude one, an exact DFT amplitude at arbitrary Fourier 
frequency r becomes 

A, = ^e-2^''-(*^-*")/^, (26) 

where Uph is the number of photons, tj is the TOA of the j^^ photon, to is the time of the start of the 
observation, and T is the total duration of the observation. Very quick harmonic-summing searches of 
an observation are possible using this technique, with the added benefit that "scalloping" (see §3.2.1) is 
non-existent. 

Since eqn. 26 only involves a summation over the number of photons, it can be computed quickly if Uph 
is relatively small. Great increases in computation speed can be had if we search a regular grid of Fourier 
frequencies. Trigonometric recurrences such as 

coa{e + S) = cos(6') - [a cos(6l) + /3 sin(6l)] (27a) 
sin(6l + 5) = sin(6l) - \a sm{e) - (3 cos{e)] , (27b) 

where a = 2sin'^ ('^/2) and /? = sin(i5), allow extremely efficient calculation of the complex exponential for 
each TOA (Press et al. 1992). This technique allows one to calculate billions of Fourier frequencies from a 
few hundred photons using only modest computational resources. 
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4. Improving the DFT Response to Arbitrary Signals 

4.1. Fourier Interpolation 

The potential for up to a 30% decrease in signal-to-noisc {S/N oc VP) due to an essentially arbitrary 
difference between the signal frequency and the integer frequency of the nearest Fourier bin is clearly a 
drawback in the use of the DFT (see §3.2.1). However, if we could calculate the Fourier Transform (FT) at 
a higher frequency resolution than the 1/T spacing that results from the FFT, we could significantly reduce 
or eliminate scalloping and effectively flatten the Fourier response as a function of frequency. 

One possibility for increasing the frequency resolution is to simply calculate the DFT by a brute force 
summation at frequencies between the integer frequencies. Such a technique is possible in special situations 
(see §3.4.2), but for most applications, the computational costs would be unacceptably high. Another well- 
known possibility is to "pad" the end of the time series with a large number of points with values equivalent 
to the mean of the data^. The padding adds no power to the data but it does increase the Fourier resolution 
since T has been artificially increased by the padding. While this technique is simple and effective for short 
time series, the difficulties involved in performing very long FFTs (§2.2) makes this technique difficult when 
dealing with long time series. 

Yet another way to calculate a higher-resolution Fourier response is to use the complex amplitudes 
produced by the standard FFT to interpolate responses at non-integer frequencies — a process known 
as "fine-binning" or "Fourier interpolation" (e.g. Middleditch et al. 1993). Similar techniques allow the full 
recovery of a signal's theoretical coherent response provided that the signal's behavior during the observation 
is either known or can be guessed. 

The purpose of Fourier interpolation is to calculate a complex Fourier amplitude at an arbitrary fre- 
quency fr = r/T, where r is any real number, such that the result is sufficiently close to the exact calculation, 

Ar=Y, n.jc-^'^'i''/^ . (28) 

We can rewrite this expression as 

N-l 

Ar=Y^Ak e-''^('-*)sinc [7r(r - k)] . (29) 

fe=0 

where the are the complex FFT amplitudes at the integer frequencies I (see Appendix A and §4.2.2 for 
a derivation and discussion of this result). 

The sine function in eqn. 29, provides the key to computing an accurate interpolated amplitude using 
a relatively small number of operations. Since sinc[7r(r — A;)] — > as 7r(r — fc) — > ±oo, the expansion of 
Ar in terms of the Ak is dominated by the local Fourier amplitudes (i.e. where k r). We can therefore 
approximate A^. as 

[r]+m/2 

~ Ak e-''^('-'=)sinc [7r(r - k)] , (30) 

k=[r]-m/2 



"Padding with the data mean is preferable to zero-padding since zero- padding introduces low frequency power into the 
Fourier response. 
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where [r] is the nearest integer to r, and m is the number of neighboring FFT bins used in the interpolation. 
Note that the interpolation is simply a correlation of the local FFT spectrum around the desired frequency 
element with a "template" response — in this case the theoretical response of a DFT to a sinusoid as 
described by eqn. 18c. 

The upper panel in Fig. 6 shows both the raw FFT power spectrum (denoted by grey dots) and the 
interpolated power spectrum (the line connecting the dots) for a radio observation of the short-period binary 
pulsar PSR J1807— 2459. The spectra cover a narrow frequency range near the pulsar's rotational frequency 
and were calculated using m = 32 and a frequency step size of Ar — 1/16 (compared to the raw FFT 
frequency step size of Ar = 1). Note that the interpolated spectrum reveals much more information regarding 
the shape of the frequency profile — including the true amplitude and location of the maximum in Fourier 
power. We will see in §5 how this can be used to deduce further information regarding the signal properties. 

A computationally less expensive version of Fourier interpolation is "intcrbinning," where we approxi- 
mate the FT response at half-integer frequencies using only the nearest 2 integer frequency bins. By using 
the Fourier interpolation equation (eqn. 30) with m = 2, ignoring an overall phase shift, and boosting the 
response such that the best-case response (at half-integer frequencies) is equivalent to the full response, we 
get 

A,^ic^^{Ak~Ak+i). (31) 

This particular formulation of interbinning was reported by van der Klis (1989), and its response is shown 
in Fig. 2. Middleditch et al. (1993) have contributed a correction to this formula for use when the data are 
padded at the end. 

Interbinning is extremely useful since such a computationally inexpensive calculation reduces the max- 
imum loss of signal-to-noise from 1 — 2/7r or ^ 36% at a frequency offset of ^ bin to ^ 7.4% at an offset of 
± (1 — 7r/4) bins. This large but cheaply-obtained reduction in scalloping can be extremely beneficial when 
searching large numbers of FFT bins and interbins. 

It is important to note that interbins as defined above have three different properties than integer FFT 
bins. First, they have different noise properties, which makes calculation of the significance of interbin powers 

much more difficiilt. Second, each interbin is correlated with the integer bins it was created from, meaning 
that interbins are not independent Fourier trials (see §2.1 for a discussion of the IFS). And finally, interbins 
do not recover the correct phase of a sinusoid at the interbin frequency. In general, since interbins are most 
commonly used during searches to simply identify signals in the power spectrum that would otherwise have 
been lost due to scalloping, these weaknesses do not degrade the usefulness of their calculation. When a 
signal is identified, a full-scale interpolation of the Fourier amplitudes around the signal using eqn. 30 allows 
accurate estimates to be made of the signal's significance and other properties (see §5). 

4.2. General DFT Response Correction 

Fourier interpolation serves as a specific example of a much more general technique — the ability to 
completely recover the fully coherent response for virtually any signal. For Fourier interpolation, we can 
exactly calculate the response of any Fourier frequency based purely on the properties of the FT. To correct 
the Fourier transform's response to a particular signal, we must know not only the properties of the FT, but 
the properties of the signal we are looking for as well. For the cases we will discuss, this ability comes in one 
of two forms: matched filtering in the Fourier domain using only the "local" Fourier amplitudes near the 
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Fourier frequencies of interest (which we call the "correlation technique"), or the straightening of the curved 
Fourier vector addition in the complex plane (which we will call "vector bending"). 

4-2.1. Correcting for Constant Frequency Derivative 

In order to illustrate these two methods, we demonstrate how to correct for a signal whose response 
is reduced due to a constant frequency derivative / (or in Fourier frequency bins, f = fT"^)- The DFT 
operates, as we noted in §3.2.1, by "derotating" the vector addition of the data in the complex plane by 
changing the phases of each of the vector elements — causing a straight line to form for a sinusoidal signal 
with integral frequency. In the presence of a frequency derivative however, the signal frequency may change 
by one or more frequency bins over the course of the observation. The complex phase corrections provided 
by the DFT will fail to completely derotate the data, and pulsation power will be "smeared" across several 
nearby frequency bins — causing a decrease in the measured DFT response (e.g. Johnston & Kulkarni 1991). 
Fig. 1 illustrates this effect in the complex plane. 

As with a frequency error, an uncorrected frequency derivative causes the vector addition to form an arc, 
although in this case quasi-parabolic rather than circular. The decreased DFT response equals the distance 
from the origin to the end of the arc. This distance is significantly shorter than that of a coherently detected 
signal which equals the the length along the arc. 

Signals with constant or nearly constant / are quite common in pulsar astronomy — especially when 
dealing with time series of very long duration. In such long observations, even the very small spin-downs 
typical of pulsars can cause a signal to drift across numerous Fourier bins. The Doppler effects of binary 
pulsar orbits cause similar frequency drifting when the observation time is much shorter than the orbital 
period. 

The "standard" method to correct for a constant frequency derivative is to "stretch" the original time 
series to compensate for the known or assumed /. This process involves re-sampling the data ensemble nj 
using a transformation similar to 

t' = t + (32) 

Jo 

where t is the time used when sampling the original data, / is the frequency derivative, and fo is the initial 
frequency of the signal. Additional details and variations on the theme can be found in Middleditch & 
Kristian (1984); Anderson et al. (1990); Johnston & Kulkarni (1991) and Wood et al. (1991). 

By stretching the data using the appropriate transform and then FFTing the corrected time series, we 
can recover the fully coherent response. Such techniques have been used with significant success in searches of 
relatively short time-series (e.g. Camilo et al. 2000). However, this technique runs into significant difficulties 
when trying large numbers of transformations using long time series where computation of the FFT is non- 
trivial. Both techniques that we will mention allow full corrections to be made to a signal without requiring 
multiple FFTs of the full original data set. 

4.2.2. Correlation Technique 

The correlation technique is the more powerful of the two methods, and uses matched filtering in the 
Fourier domain to "sweep-up" signal spread over a number of frequency bins into a single bin. In astrophysical 
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applications, we usually have some sort of "pure" signal (like a harmonic from a millisecond pulsar), whose 
frequency changes as a function of time due to some other process (such as orbital motion or pulsar spin- 
down). In the Fourier domain, these processes cause the perfect sine- like response of a harmonic to be spread 
over numerous local Fourier bins — in effect, the sine response is convolved with a finite impulse response 
(FIR) filter (where finite in this case refers to a small portion, say m bins, of the frequency range analyzed 
rather than a short period of time). If we can predict the complex form (and phase) of that FIR filter, we 
can recover the coherent response (i.e. the perfect sinc-function) by correlating the appropriate Fourier bins 
with a "frequency-reversed" and complex- conjugated template that matches the filter. 

In mathematical terms, consider a signal with a normalized Fourier response of Ak-ra^ where k — Vo 
is simply the frequency offset of bin k from some reference frequency To, which goes to zero as |fc — ro| 
approaches some number of bins m/2. For Fourier interpolation as described in §4.1, this response is equal 
to eqn. 18c without the Aq factor (i.e. normalized to an amplitude of one for a coherent response). The 
complex-valued Fourier response of such a signal at frequency To can be calculated with the sum 

[ro]+m/2 

Yl ^>o^l-k- (33) 

fe— [ro] —rn/2 

If To is initially unknown (i.e. we are searching for a signal with the response shape as defined by the template 
but at an unknown frequency) we simply compute this summation at a range of frequencies r. 

Calculating eqn. 33 over a range of evenly spaced frequencies is equivalent to correlating the raw FFT 

amplitudes with the template and is therefore most efficiently computed using short FFTs and the convolution 
theorem. With FFTs of length Af, such that m <C M <C N/2, we can search a very-long FFT of length N/2 
for any signal whose Ak-r^ we can compute, using overlap-and-save or overlap- and- add techniques (see e.g. 
Press et al. 1992). Such calculations have advantages over standard time-domain "stretching" techniques in 
that that they are memory local and can be easily parallelized important properties when dealing with 
very-long time-series and modern distributed memory computer architectures. 

Moving to our example of a signal with a constant frequency derivative, a single harmonic of the signal 
has the form 



n{u) = a cos 



27r ( roU +-u']+> 



(34a) 
(34b) 



_ ^ g2-n-i(roM-F§u^)gi(/> _|_ g-27ri(r„u-|-§u^)g-i 

where we use the same notation as §2.1. Neglecting the second term as in §3.2.1 and Fourier transforming 
at some "center" or average frequency = r -|- r/2, we get 



(35) 



where qr = Tc — r^, and the real "center" frequency of the signal is Vc = To + r/2. This integral can be 
evaluated in closed form 



L 



1 12 

gi^(r«^+2g,«) ^ e-i^^ {5 ^Zr) - S (Yr) + i [C (Yr) - C (Zr)]} (36) 



/o V2r 
where Yr — \j2/fqr, Z^ — \/2/f [q^ + r) , and C{x) and S{x) are the Fresnel integrals 



C{x) = cos (1*2 j dt, S{x) = £ sin (^t^ j dt. (37) 



The Fourier transform response then becomes 




{S {Zr) - S {Yr) + i [C {Yr) - C {Zr)]} . 



(38) 
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Using the correlation technique, the coherent response can be recovered by convolving local Fourier 
amplitudes with the "frequency-reversed" and complex-conjugated template as defined by eqn. 38. This 
response, at average Fourier frequency Tc and Fourier frequency derivative r can therefore be written as 



Eqn. 39 takes into account the fact that the signal has been "spread" relatively evenly into the f closest 
frequency bins to Tc, while an additional small amount of signal has "leaked" into bins further away — much 
like the non-zero wings of the sinc-like response to a constant frequency signal. As a rule of thumb, the 
correct Fourier amplitude will be well-approximated if m is chosen such that f < m < 2f. 

Large scale searches of pulsations with constant frequency derivatives have been conducted using the 

correlation technique. A successful example is a search for pulsars in globular cluster NGC 6544 using 
data taken with the Parkes radio telescope (Ransom et al. 2001). The search was conducted using an FFT 
of 13865600 points over Fourier r values from —100 to 100 with a step-size of Af = 2 and included the 
calculation of amplitudes at half-bin frequency intervals. The search resulted in the detection of the 3.06 ms 
PSR J1807-2459 in a low-mass binary with orbital period 1.7 h, the second shortest radio pulsar orbital 
period known. A detailed view of the pulsar's fundamental harmonic is shown in Fig. 6. The plot was 
calculated using the correlation technique with spacings of Ar = 0.0625 and Af = 0.25. The generation of 
such a piece of the /-/ plane takes only a fraction of a second on a rather modest workstation. 



Vector bending is one of the simplest and most straightforward methods to correct a Fourier response 
that has been smeared over several local frequency bins. As we described in §2.1, the DFT can be thought of 
as the vector addition of N complex numbers. This addition produces a straight line in the complex plane for 
a coherently detected sinusoid. For a sinusoid with a non-integral or time-varying pulsation frequency, the 
standard DFT addition produces a curved shape (see Fig. 1). Since the amplitude of the Fourier response 
is the distance between the origin and the end-point of the vector addition, any curvature in the vector 
addition implies non-optimal signal detection. 

The precise shape of the response curve in the complex plane depends on the mismatch of the signal's 
(possibly time-dependent) pulsation frequency and the frequency used in forming the DFT addition (i.e. the 
closest FFT bin). Regardless of the shape, though, for short enough segments of the curve, the segments 
differ little from straight lines. We can therefore approximate the shape of the curve as a sum of G linear 
segments, each of which contains N/G points from the full-resolution vector addition. In terms of the r**^ 
DFT amplitude, we can write this as 




k=[r]-m/2 



(39) 



4-2.3. Vector Bending 




(40) 



9=0 



g=0 h=0 
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where j = gN/G + h. This is equivalent to calculating and then summing G independent DFTs (the Br,g), 
each of which suffers virtually no loss in sensitivity when the curvature over a segment is small. 

Using these vector addition segments or sub-vectors, we can correct for the loss of sensitivity due to 
curvature by simply straightening the vector addition. If we can predict the true pulsation frequency of a 
signal as a function of time, we can predict how much curvature will accumulate in each sub-vector and then 
remove it by rotating the segment appropriately. 

For our example of a signal with constant frequency derivative, the instantaneous phase (ignoring the 
intrinsic phase of the signal) is equal to 

^true{u) = 2tt j r{u) du (41a) 

= 2tt j {vo + fu) du (41b) 
= 2-KroU + irfu^, (41c) 

where rg = foT is the initial pulsation frequency of the signal, and r = /T^ is the frequency derivative. The 
process of taking a DFT removes an instantaneous phase equivalent to ^dft{u) = 2-Kru from the signal (see 
eqn. 2). So the instantaneous phase error is equal to 

terror («) = ^true{u) - ^DFt{u) (42a) 

= 2'K{ro - r)u + wru^ . (42b) 

Therefore, to correct a particular signal using vector bending, we first calculate the Br,g using eqn. 40 for 

a particular Fourier frequency r (such as the frequency of a known pulsar). Now we attempt to un-bend 
the full Fourier response by summing the B^-^g after correcting for the phase errors ^ error {ug) as defined by 
eqn. 42a. The corrected response is equal to 

G-i 

Ar,f = S^.se-'*''-'-'"-^"^^ (43) 

3=0 

A choice of G ^ 10^ will essentially eliminate the loss of response for reasonable frequency offsets and 
frequency derivatives (i.e. less than a few 10s of Fourier bins). 

While impractical for large-scale searches due to the fact that the _B,..g must be recomputed every few 
r, vector bending offers significant computational advantages in certain situations. In particular, X-ray 
observations often consist of short (< 1 hour) "on-source" segments separated by hours, days, or even weeks 
of "off-source" time (see also §3.4). An FFT of the entire time series might be prohibitively expensive. 
However, if we can determine an "initial guess" frequency (e.g. by FFTing one segment of the observation 
or from an ephemeris), we can quickly calculate the Bk^g at this frequency from the "on-source" intervals 
alone. We can reconstruct the /-/ plane around our frequency of interest without ever creating the full 
"fiUed-in" or padded time scries, let alone calculating a potentially huge FFT. Such techniques have allowed 
us to perform frequency analysis of very-long stretches of data from ROSAT observations of PSR B0540-59 
(Eikenberry et al. 1998). Similarly, Fig. 1 shows the use of the method for a 2.4 day observation of the Crab 
pulsar. 
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5. Signal Property Estimation 

Besides correcting for losses of sensitivity, Fourier interpolation and the other response correcting tech- 
niques mentioned in §4.2 allow us to determine other useful properties of a detected signal. Using detailed 
amplitude and phase information from a signal's Fourier harmonics, we can estimate properties such as the 
statistical significance of the signal, the location and duration of the signal in the time series (the "centroid" 
and "purity" respectively), the precise pulsation frequency, and estimates of the measurement errors for 
Fourier power and phase. 

The first step when estimating signal properties in the Fourier domain is to isolate the true peak of 
the Fourier response in power. This is easily accomplished by using the matched filtering techniques to 
generate an oversampled grid of amplitudes near and around the signal candidate (see Figs. 7 and 6). Simple 
optimization algorithms such as the downhill simplex method can then be used to refine the peak location 
(e.g. Press et al. 1992). Once the peak has been located, estimates of the first and second derivatives of 
power and phase with respect to Fourier frequency, obtained using Fourier interpolation, can be used to 
calculate various useful signal properties (Middleditch et al. 1993). 



5.1. Power, Phase and Signal Amplitude 

When the peak of the Fourier response has been located as a function of Fourier frequency and the 
other search parameters, the measured power is defined as 

\Ar f 

Pmeas ~ j ('^4) 

^norjn 

where Pnorm is the expected noise power and is usually described by one of {d^)^ Piocah or Hph as discussed 
in §3.1. Groth (1975, see §3.3) showed that since the measured power is a random variable due to the presence 
of noise, its variance is ^Pgignai + 1, where Psignai is the power caused by the signal. Since we do not a priori 
know the true signal power, a good estimate for the variance of the measured power is simply 2Pmeas — 1 

since {Pmeas} Psignai 1- 

Using Pmeas 1 as well as the knowledge that a sinusoid of amplitude a in a noisy time series produces a 
power with an expectation value of {Pmeas) = a'^N^/ {^Pnorm) + 1 (see §3.2.1), we can estimate the signal 
amplitude as 

2 ; 

{a) = J^V Pnorm {Pmeas - !)• (45) 

For binned data containing a signal with Fourier frequency r, the measured power should be multiplied by 

l/sinc^(7rr/A^) to correct for the loss in sensitivity due to binning (see §3.4.1). Vaughan et al. (1994) provide 
detailed instructions on how to estimate upper limits on pulsation amplitudes as well as estimates of the 
overall sensitivity of a search. 

The statistical significance of a signal is also determined by Pmeas ■ The probability that noise caused a 
particular power to equal or exceed Pmeas is given by e"^"**"*^ (eqn. 10 with P' = Pmeas)- But for a search 
over Njps independent Fourier powers, the probability that at least one of the noise powers exceeds Pmeas 
is given by 

Prob{Pnoise > Pmeas) = 1 - (l - e"^--^ )'^'^^ . (46) 

Vaughan et al. (1994) show how to use this information to set detection thresholds that minimize the number 
of spurious candidate signals and give high confidence that signals with powers above the thresholds are real. 
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Using the real and imaginary parts of the peak Fourier response, we can also calculate the phase of the 
sinusoidal signal as 

llm{Ar,...y 
_Re{Ar,...)_ 

radians. Using similar arguments as for the measured power, the variance of the measured phase is approx- 
imately l/{2Pmeas — 1) radians. 



= arctan 



(47) 



5.2. Signal Location and Duration in Time 

Astronomical observations of pulsations effectively consist of a window of on-source time where pulsa- 
tions are present, and the rest of the time when they are not. For most of this paper we have assumed that 
a signal is present throughout the observation as evidenced by the limits of integration for cqn. 2, which in 
time-normalized units go from < u < 1, or equivalently, from <t <T. In effect, pulsations such as that 
defined in eqn. 14a are multiplied by a square window function defined as 1 during the observation and at 
all other times. This window function is simply a property of the DFT and is due to the finite duration of 
our observation. 

It is possible, though, and for various reasons often likely, that a signal we are observing turns on and 
off or varies in intensity during an observation. The behavior of the signal itself effectively defines a new 
window function, W{u). By measuring the moments of this window function with respect to time, we can 
determine approximately where in our data a signal is located and for how long. 

The approximate location of a signal in a time-normalized time scries is described by its centroid, 
C — (u) — {t) /T which is proportional to the first moment of the window function with respect to time. 
More specifically, 

g= ;s • (48) 

/_oo ^(^) 

Middleditch & Cordova (1982) wrote this in terms of the measured Fourier response as 

27r dr ^/24P{ro) 

where P{ro) and ^{ro) is the phase measured at the peak of the Fourier response (i.e. (pmeas) and r is 
the Fourier frequency (see Appendix D for a derivation). Signals present throughout an observation have 
(7 = 1/2 while those present in only the first or second halves of the observation have (7 = 1/4 or (7 = 3/4 
respectively. 

The second moment of the window function with respect to time is related to the moment of inertia of a 
function and can therefore be used to estimate the root-mean-squared (rms) dispersion of the pulsations in 
time about the centroid. Using this information, Middleditch & Cordova (1982) defined a parameter called 
the "purity" (and symbolized by a) as 




/ , (50) 



2P{ro) dr^ a^lOP{ro) 

where P{ro) is the measured power at the the peak of the Fourier response (i.e. Pmeas, see Appendix E for 
a derivation). The scaling in eqn. 50 is chosen such that the rms dispersion of the signal about the centroid 
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for a window function W{u). is equivalent to that of a rectangular window function of duration a (in units 
of the time series length) centered on the centroid. A signal present throughout the data would have a = 1, 
while one present in only half the data (in a continuous section) would have a = 1/2. Signals present only 
at the start and end of an observation but absent in the middle have a > 1. Purity can also help to identify 
sidelobes caused by a periodic modulation of a signal as these Fourier amplitudes have a = \/3. 

Since the location and duration of a signal in an time series affects the Fourier response, it is important 
to understand how eqn. 18c changes if a signal is present during only part of an observation. In Appendix F, 
we show that when close to the peak of a signal's Fourier response, 

Ar ~ Ao e-2'^'C'('-''°) sine [7ra(r - Vo)] (51) 

where Ao = Na/2 c'*°, (j)o is the intrinsic phase of the signal, and Tq is the true signal frequency (in FFT bins). 
This equation demonstrates that for centroids different from 1/2, the phase shift between consecutive FFT 
bins differs from the n radians shown in eqn. 18c. Similarly, for purity values different from 1, neighboring 
FFT bins show more or less correlation with each other (i.e. the central peak of the sine function changes 
its width). 



5.3. The Pulsation Frequency and Frequency Derivative 

The true pulsation frequency of the signal is located at the point where dP/ dr = 0. Furthermore, given 
the response in eqn. 51, we can show (Appendix B) that the uncertainty in this measurement (in Fourier 
bins) is given by 

= 7^=- (52) 

This uncertainty is considerably smaller than the often-quoted "frequency error" for the FFT of one bin 
width, which is simply the frequency resolution returned by the FFT algorithm. 

If the correlation method is used to isolate a peak in the /-/ plane as shown in Figs 7 and 6, we 
can calculate the uncertainty in the measured r value by using similar arguments and methods as for the 
frequency uncertainty (see Appendix C for a derivation). The uncertainty in the r (in Fourier bins) is 
approximately 

3^10 (53^ 



6. Conclusions 

In this paper we have described techniques that allow sophisticated and fully coherent Fourier analysis 
of very long time series. Most of these techniques use the wealth of information provided by the Fourier 
phases — information discarded during "standard" analyses based on raw power spectra. 

Significant gains in sensitivity and efficiency are possible when using Fourier phase information during 
the search for periodic signals (using the Fourier domain matched ffitering techniques described in §4.2) and 
when characterizing signals that are known to be present in the data (using the parameters described in 
§5). The methods of Fourier domain matched filtering allow efficient, memory local, and inherently parallel 
analysis of extremely long time series with only modest computational resources. Billion point out-of-core 
FFTs followed by fully coherent matched ffitering pulsation searches are possible on "standard" workstations. 
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More traditional time domain based techniques (such as acceleration searches performed by stretching or 
compressing the time series followed by large in-core FFTs) on similarly sized time series require specialized 
high-performance computing resources, assuming they can be performed at all. 

As astronomical instruments become more sophisticated and specialized, time series of ever increasing 
duration and time resolution will appear. The Fourier domain techniques described in this paper should 
prove to be essential tools in their analysis. 
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A. Derivation of Fourier Interpolation 

Following the derivation found in Middleditch et al. (1993), we begin with the definition of the k*'^ DFT 
element 

Xjn.e-^-'^'^/^, (Al) 

3=0 

which we then rewrite by substituting the inverse DFT for the nj 

j=0 \ 1=0 J 



1=0 j=o 

The last summation can be computed exactly using the identity 
such that when N ^ 1 we have 

JV-l 

^-2nij{k-l)/N ^ ^-i^(fe-()(l-i: 

,.=0 sm[wik-l)/N] 

^ i^(k-i) sin Hk - I)] 
T:{k-l)/N 

^ iVe-'''('=-')sinc[7r(/e-0]- (A4c) 



V- -2nij{k-l)/N _ i,r(fe-0(l-i) jMiiiLlllL ( A4a] 

(A4b) 



Substituting this expression into eqn. A2b and changing the integer frequency k into a continuous real-valued 
frequency r, we arrive at eqn. 29 



JV-l 



Ar=J2'^i e-''^(''-')sinc [7r(r - /)] • (A5) 



(=0 
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B. Derivation of Frequency Uncertainty 



For a given Fourier frequency offset A^. = r — ro, where Vo is the Fourier frequency of the signal, the 
magnitude of the Fourier response and the power go as (see Appendix F) 



\A{r)\ = |Ao|sinc(7raAr) 
P{r) = Pmeas sinc^ (naAr) 



(Bla) 
(Bib) 



where a is the signal purity. We can expand the sine function in order to approximate the expression for 
the power near the peak of the response as 



P{r) = Prr 



sin(7raAr) 


2 

^ P 


TraAr 


— ^ meas 



1 - 



(TraAj.) 



Taking the derivative of power with respect to r and solving for A^ we obtain 

3 dP 

A, = - 



dr 



As expected, when the Fourier frequency equals the true frequency of the pulsations (i.e. A^ 
response peaks and dP{r)/dr = 0. 



(B2) 
(B3) 

0), the Fourier 



In order to estimate the uncertainty in A^ , we apply standard propagation of errors to arrive at 

3 



27r2a2R, 



-crp/(r). 



(B4) 



where we have replaced dP{r)/dr with P'{r) to simpliiy the notation. The derivative of the power at the 
true frequency can be approximated using finite differences as 



P'(r) 



P{ro + K)-P{ro-K) P+-P- 



2Ar 



2Ar 



(B5) 



where we have simply renamed P {vo + A^) and P {tq — A^). The uncertainty in P'{r) can also be approxi- 
mated using finite differences and error propagation. Since P"*" and P~ are highly correlated when A^ -C 1, 
their uncertainties are also correlated giving 



''^'('-)^2a:^''^"+''^-^^ a7 



(B6) 



Now, we turn to the question of the uncertainty in P+, closely following Middleditch (1976). The 
amplitude of the Fourier response at the true signal frequency can be written as 

V Pmeas = Vj COs{^j) (B7) 

j=0 

where yj = nj/\/Pnorm, are the points in our time series as defined in eqn. 1, but scaled using the appropriate 
Pnorrn such that the measured power, Pmeas, is properly normalized (see §3.1). Similarly, the (pj represent 
the pulsation phases at times u = j/N, but rotated by the measured Fourier phase, <pmeas, such that the 
result of the vector addition lies along the real axis in the complex plane (i.e. the final complex phase is 
zero). In effect, this transform isolates components of the data that are parallel to the final Fourier response. 



- 23 - 



At a small frequency offset from the true frequency, we can expand the power in a similar fashion as 

N-l 

Vp+= ^y,.cos(0, + 5^^), (B8) 
i=o 

where 6^^ are the "phase errors" introduced by the frequency offset. The phase errors add curvature to the 
vector addition and are defined as 

6^^ = 27raA^ (j- - = 2naAr " ^ > (B9) 

where u is the normalized time, u = t/T = j/N. The 1/2 term in eqn. B9 removes the accumulated phase 
error over the course of the observation (i.e. 27rQ!A,.u du = TraA^) and makes the vector summation of 
eqn. B7 finish on the real axis. Expanding the cosine in eqn. B8 gives 

N-l N~l 

\/P+ = ^ yj cos{4>j)cos {5^.) - ^ Vj sin {5^-) . (BIO) 

3=0 j=o 

Considering the uncertainties in the separate terms of eqn. BIO, since the cosine term is an even function 
of 5^^ , it is symmetric about Tq, and therefore does not contribute to the uncertainty in the P'{r) measurement 
as defined by eqn. B5. For the sine term, given that the cos(^j) "derotates" the signal onto the real axis by 
definition, we see that 

N-l 

^%sin(</,,) = 0. (Bll) 
j=o 

Furthermore, since the 0j and 5^^ are uncorrelated, the average value of this term will be zero, 

^ E ) sin (^^, ) ^ = 0, (B12) 

and has no systematic effect on P+. However, we can calculate the fluctuations introduced by this term 

Ew sin (</.,) sin (5^.) j ^ = (j^yjsm{cl>^)^ ^ (sin^ (5^ . ) ) . (B13) 

where the cross-terms average to zero since (pj and 5^. are uncorrelated. Due to the normalization of the j/j, 
the sum component averages to 1/2. The 5^^ component has an average of 

(sin2(5^J) = [ sin' {6^^) du (B14a) 

Si du (B14b) 







Jo 





1 

2 I „,2 



1 ^ 1 



(2TTaArY l^u' -u+ -j du (B14c) 

(B14d) 



-24- 



Therefore eqn. B13 is equal to (7raAr)^/6, and the variance of ^/P+ will be 



a\ = ^JP+ 



6~~ 



(B15) 



To get the standard deviation we take the square root of this expression. Propagating errors to get the 
uncertainty in P+ adds a factor of 2^/P+ to give 



p(r o+£^r) = 27raAr 



Prr 



(B16) 



where we have used the approximation P+ ~ Pmeas 
have 



. Finally, substituting into equations B6 and B4, we 
^ (B17) 



^Ta^/6P^^ 

A much simpler derivation of eqn. B17 is possible if we realize that properly normalized powers times two 
are distributed according to a distribution with 2 degrees of freedom (see §3.1). For such a distribution, 
a la error corresponds to a change in the measured of 1/2. In the case of powers, the zLla errors can 
therefore be found by starting with the expansion of the power around P{ro) = Pmeas as given in eqn. B2, 



and then solving for the A^ that corresponds to P„ 



,os - 1/2. This gives 
3 

Tray/6Pmeas 



(B18) 



(B19) 



C. Derivation of / Uncertainty 

If a peak in the / / plane has been isolated using techniques similar to those shown in §4.2.2, we can 
calculate the error in the measurement of the true Fourier frequency derivative r'o = foT^ in a manner 
similar to that for the frequency uncertainty as described in Appendix B. Signals with non-zero frequency 
derivatives have Fourier peaks which are located off the r = line in the / / plane, but the shapes of those 
peaks are independent of f and in fact depend only on the window function of the signal (see Appendix F). 
The shape of the response in power as a function of Aj. = r — r'o at the "correct" Fourier frequency r = To 
is described by eqn. 38, and can be written as 

when Qr is defined as 

Qr = r,~r^=[ro + Y)-['o + ^ J = 

This definition of Qr keeps the magnitude of the Fourier response symmetric about 
value of A^. The power as a function of A^. is therefore 




(C2) 

Qr no matter what the 
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If we expand the Fresnel integrals about as 



C{x) ~x- . . . , and S{x) ~ . . . , (C4) 



and then substitute, the power becomes 



1 - 



(C5) 



(C6) 



180 

Taking the derivative of P{ro, r) with respect to f and solving for Aj. gives 

90 dP{ro,f) 
■K^a^Po df ' 

Prom propagation of errors, 

90 

where P'{ro,r) = dP{ro,r)/dr. Similarly, after using a finite difference estimate of the power derivative 
following Appendix B, we find 

^P'{ro,r) — T 1 {^°) 

where P{ro, A^.) represents the power as measured at f = r'o + Aj.. 



Closely following Appendix B, we represent sJPiro-, Af) as a sum of the parallel components of the 
properly normalized time series points, y^, as 

^P{ro,^f) = ^ Vi cos + ) . (C9) 

i=o 

The 8^. are the "phase errors" introduced when A^ 7^ 0, and are defined as 

5^^. = TTQ^A^ ( - M + M , (CIO) 



where the u term comes from keeping the response symmetric about (i.e. g,. = —o? /2) and the 
1/6 removes the accumulated phase error over the course of the observation (i.e. /g Tra^Aj. (w^ — u) du = 
— 7ra^Aj./6) and makes the vector summation of eqn. C9 finish on the real axis. 

Expanding the cosine term of eqn. C9 gives 

N-l N-1 

^/P{ro,Af) = ^ Vj cos (0j) cos (J^J - ^ yj sin (^j) sin (6^^) . (Cll) 

The first term shortens both P{ro. A,- ) and P{ro, — Aj.) by the same amount and docs not affect the derivative 
of power. For the sine term, since and are uncorrelated, the average value is zero (see Appendix B), 
but its fluctuations are important. To calculate the fluctuations, we square the terms and get 

^^2/,sin(0,)sin(<5^jj ^ ^ UY^Jj^^^i^U {sin^S^^)) . (C12) 
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Due to the normalization of the yj, ^ ^X^jLo^ Vi (^i)) j averages to 1/2 as before, and we can directly 
calculate (sin^ i^4>3)) ^ 

(sin2(^0,)) = / sin2(5^.)dw (C13a) 
Jo 

- [ S^du (C13b) 
Jo 

= n'^a^Aj £ (u'^-u+ du (C13c) 

= (C13d) 

The fluctuations from the sine term are therefore (7ra^Ar)^/360, and since squaring eqn. C9 doubles the 
errors, the standard deviation of P(ro, A^) is 

27rQ;^Ar f-= na^Af 



(^P(r„,A.) = VP{ro,Af) —f=^ ^ \J Pmeas (C14) 

VoOL) v90 

Substituting into equations C8 and then C7 as in Appendix B gives us the uncertainty in the frequency 
derivative 

" " V meas 



D. Derivation of Centroid 

The centroid is a measure of the approximate location of a signal in a time series as estimated by the 
first moment of the signal with respect to time (see §5.2). We can think of a sinusoidal signal in our data 
as being always present but modulated in intensity by some window function W{u), wlic;rc u = t/T is the 
normalized time and T is the length of the observation. A "normal" observation of a pulsar of constant 
intensity would therefore have W{u) = 1 when < u < 1 and W{u) = at all other times (i.e. a square 
window). Our signal is therefore described by 

s{u) = a cos {2'!TroU + 4) o)W{u), (Dl) 

where a is the amplitude, Tq = foT is the Fourier frequency, and (j)o is the phase of the sinusoid at time 

u = (). 

Since the centroid of a function is proportional to the first moment of the function with respect to time, 
we can easily calculate the centroid using the Moment Theorem of Fourier transforms. Bracewell (1999) 
does this and defines the centroid as 

where A{Qi) and ^'(0) are the Fourier transform and its first derivative with respect to r measured at r = 0. 

Eqn. D2 is not directly applicable for our sinusoidal signal since the information about the window 
function in eqn. Dl has been shifted to the frequency of the sinusoid in accordance with the Modulation 
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Theorem of Fourier transforms. Accordingly, we can apply the Modulation Theorem to eqn. D2 which gives 
us 

2niA{ro) 



Finally, we can write A{r) in phasor form as A[r) = a{r)&'^^^^ where a(r) and (f>{r) represent the Fourier 
amplitude and phase as functions of the Fourier frequency r. The derivative with respect to Fourier frequency 
can be written 

A'ir) = ^e'^(^) + a{r)i^e^^^^\ (D4) 
or or 

At the frequency of our signal, the amplitude is A(ro) = a(ro)e'"^('"°) and the Fourier response is at its peak, 
making da(ro)/dr = 0. Therefore 

A'{r„)^A{r,)\^^. (D5) 
Substituting eqn. D5 into eqn. D3 we arrive at 

(D6) 

which is equivalent to eqn. 49. 

We can also estimate the uncertainty on the measured value of the centroid. Following Appendix B from 
eqns. B8— B15, we note that the same noise fluctuations introduced when offsetting from the true pulsation 
frequency by a small amount that effect the Fourier amplitude will also effect the Fourier phases. When 
Ar is small and the powers are properly normalized, the amplitude fluctuations of eqn. B13 with variance 
7r^A^/6 correspond to an uncertainty in the phase measurement of cr0(r^+Ar) — 7rAr/-\/ 6P(ro) radians and 
therefore 

1 . (D7) 



E. Derivation of Purity 

The "purity" of a signal (see §5.2) is a measure of the rms dispersion of the pulsations in time with 
respect to the centroid and is directly proportional to the variance in time of the window function of the 
sinusoid from eqn. Dl. The time variance is defined as 

but can be written as 



{u-{u)y) = {u^}-{uy (E2a) 



A"(0) , [A'm 



2 



47rM(0) + 47r2[A(0)]^' ^^^^^ 

using the Moment Theorem for Fourier transforms (see e.g. Bracewell 1999). Since our signal is sinusoidal 
(see Appendix D), application of the Modulation Theorem gives 

(u-{u)) ) = — ^ „], — rH — ^ — (E3a) 
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Using A{r) = a(r)e"^('"^ and remembering that da{ro)/dr = 0, the second derivative of the Fourier 
amplitude is 



dr'^ \ dr J a{ro) dr'^ 

Prom Appendix D, we know that d(j){ro)/dr = —2itC which makes the first term of eqn. E4 equal to zero 
since d'^^{ro)/dr'^ = and the second term equal to —A'K'^A{ro)C'^- The second derivative of power at the 
peak response can be written as 

= 1^ {A*{rM{ro)) = 2a(r„)^, (E5) 
making the third term equal to 2(a(r°)V ^ ' Substituting into eqn. E3a and simplifying gives 

/.W = - 8.^P(r„) dr^ • (^'^ 

If we normalize the variance using the value obtained for a signal present throughout the observation 
(i.e. a square window where W{u) = 1 from < m < 1 and zero elsewhere — which we will call a pure signal) 
where 

(«-M''U = <"^>-<»>^ = 5-(5)' = n- 

and then taking the square-root, we are left with 



\ \ / f 



pure 



which is equivalent to eqn. 50. 

In order to estimate the uncertainty in the measurement of a, we note that by squaring eqn. E8 and 
substituting the finite difference approximation of 



(E9) 



d^Pjro) P{ro + Ar) + P {ro - A^) - 2P (r^) _ P+ + p- - 2P (r^) 

where A^ corresponds to a small frequency offset from the measured peak power P(ro), we get 

3 P++P--2P(r,) 

This means that ^ 

''"^47r2aA2p(r„)''^- ^^^^^ 
where Ap = P+ + P~ — 2P (ro) and the extra factor of 2a comes from converting to ct^. 

We can then expand the Fourier amplitude around the peak of the signal as in eqn. BIO, where we 
see that the amplitude fluctuations from the sine term (which is anti-symmetric about Tq) will cancel after 
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the addition of P{ro + A^) and P{ro — A^) in the finite difference approximation shown above. Conversely, 
the fluctuations due to the cosine term (which is symmetric about Tq) will add. These fluctuations can 
be computed by taking the variance of the first non-constant term in the expansion of cos{5j) which is 
S]/2 = 27r2A2(u - 1/2)2. rpj^g computation gives 



('^1/2 -(^1/2))') = ((<5|/2)^)-(^|/2)^ 





{l 




n'At 


1 


20 ~ 


n^Af 




45 





u - - 1 du 



du 



(E12a) 
(E12b) 

(E12c) 
(E12d) 



and since the variations in the cos{4>j) term average to 1/2, just like the sin(0j) term did in eqn. B13, the 
standard deviation of ■^/P{ro + Ar) is equal to 



90 



90 ■ 



(E13) 



Doubling the error when converting from amplitudes to powers and then doubling it again due to the addition 
of the errors from the power offsets in eqn. E9 gives 



(jAp = A'k'^AI 
Finally, by substituting into eqn. Ell, we get 

-477^ A: 



90 



4n'^aA'^P{ro 



P{ro) 



90 aVlOP(ro) 



(E14) 



(E15) 



F. Centroid, Purity, and Fourier Response 

In order to consider the effects of centroid and purity on the Fourier response to a sinusoidal signal as 
described by eqn. Dl, we initially assmnc a Fourier response equal to eqn. 18c of 

Ar = Ao e-''^('-''<')sinc [7r(r - r^)] (Fl) 

where Ao = Na/2 e"'^", (po is the intrinsic phase of the signal, and Tq is the true signal frequency (in FFT 
bins). This response is correct only for signals with a square window function (i.e. W{u) = 1 from < m < 1 
and zero elsewhere). 

From eqn. Fl we see that a change in Fourier frequency of a single Fourier bin causes a change in the 
measured Fourier phase of tt radians. This phase change is also visible from the centroid equation for a pure 
signal with C = 1/2, where d(/)(r) = — tt for every dr = 1. Rewriting the centroid equation and integrating, 
we see that the Fourier phase near the peak response goes as 

(pir) = -2-KrG + c. (F2) 
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When r = To, (I) = (j)o, allowing us to solve for the constant of integration, c = <po + '^'kToC. Substituting into 
eqn. F2, we see that the phase of the Fourier response is equal to 

-^(ro) =(l>o- 2wC{r - r^). (F3) 

Therefore, for signals that have centroids different from 1/2, the phase change across a single Fourier bin is 
different from the usual tt radians. 

The purity, as described in §5 and Appendix E, is the cfFoctivc duration of a square window which 
reproduces the measured rms dispersion of the signal in time about the ccntroid. Since the Fourier response 
to a square window goes as sinc(7r/r), where fT = r, we can see that replacing the window of length T with 
one of effective duration aT causes the Fourier response to go as sinc(7ra/T). This fact is also approximately 
true for more complicated window functions as long as |r — rd <C 1. The Fourier response to a windowed 
sinusoid is therefore 

Ar ~ Ao c-2'^iC'(r-ro) gjj^^ [^^(^ _ (p4) 

Numerical simulations show that this approximation is valid for purity values a < 1.5. 

These same "effective duration" arguments also apply to the shape of the response in the / direction 
of the /-/ plane. A change in the effective duration of a signal causes a change in the / response since 
/T^ — > /a^T^, or equivalently, r ro? . The results of this change can be calculated by directly substituting 
fo? for r in eqn. 38, as was done in Appendix C. 
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Frequency ('r') Offsets Frequency Derivative ('i^) Offsets 




PSR J1 807-2459 Orab Pulsar with ROSAT 

Fig. 1. — Fourier responses plotted as a series of vector additions in the complex plane. The outer circles in each plot show 
the Fourier amplitude of a signal where all power is recovered by the vector addition (i.e. calculation of the DFT at the correct 

signal frequency r and frequency derivative r for signals with linear changes in frequency over time). The end-points of the 
vector additions are the Fourier amplitudes. For plots (a) and (b), a fully recovered signal would start at + Oi and end at l + Oi. 
(a) and (b) show the effects on Fourier amplitude and phase when a signal's intrinsic frequency {r in bins or wavenumber) or 
frequency derivative (r in bins/observation) differs from the computed values. For plot (b), the average Fourier frequencies in 
each case were correct, and only the frequency derivatives were in error, (c) shows the response of PSR J1807-2459 during its 
discovery observation (Ransom et al. 2001) with and without corrections for pulsar acceleration (r) and interpolation in Fourier 
frequency (r). The vectors were calculated using the method shown in §4.2.3. The fact that the corrected vector does not quite 
reach the circle implies that higher order effects of the orbital motion remain un-corrected (see Fig. 6). (d) shows corrected and 
uncorrected responses of 10,000 randomly selected photons from a 2.4 day ROSAT observation of the Crab pulsar. 
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0.2 0.4 0.6 0.8 1 

Fourier Frequency Offset from Integer Bin 

Fig. 2. — The thin solid black lines plot the sine function amplitude responses of raw (or integer) FFT bins. The point where 
these lines cross (at an offset of ^ and an amplitude of 1 — 2/7r ^ 0.637 of the full response) is the worst case response for an 
uncorrected FFT (see eqn. 18c). The thin dashed line is the response of an "interbin" as calculated using eqn. 31. The thick 
black line shows the overall "scalloping" when interbinning is used. Worst case responses with interbinning occur at offsets of 
± (1 — 7r/4) ~ 0.215 with amplitudes of ~ 0.926 of the full response. 
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Pulse Phase 



Fig. 3. — Sample pulse profiles from the modified von Mises distribution (MVMD) as described in §3.2.2. FWHM is the 
fractional full-width at half-maximum and k is the MVMD shape parameter. High values of k result in Gaussian profiles, while 
as K — » 0, the pulse shape becomes more and more sinusoidal. The integral of a full pulse is equal to one unit, all of which is 
pulsed (i.e. the pulsed fraction is one). 
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Pulsation Duty— Cycle or Fractional FWHM 

Fig. 4. — The thick soUd line represents the approximate number of harmonics from an MVMD signal (see §3.2.2 and Fig. 3) 
that produce Fourier amplitudes greater than one-half the amplitude of the fundamental. The thin solid line is the 1 /FWHM 
rule-of-thumb that is often used to estimate the number of significant harmonics a signal will generate. The thin dashed line 
plots the ratio of the fundamental amplitude for an MVMD signal to a sinusoidal amplitude of the same pulsed intensity. The 
grey histogram shows the distribution of pulse widths (FWHM) for over 600 radio pulsars. 
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Fig. 5. — The thin lines represent sensitivities to MVMD signals (see §3.2.2 and Fig. 3) for the incoherent summing of 1, 2, 
4, 8, 16, or 32 harmonics, as compared to searches made without harmonic summing (§3.1). Lower numbers represent better 
sensitivity (i.e. fainter signals are detectable). The best possible sensitivity using incoherent harmonic summing is shown by the 
thick black line. It should be noted that incoherent summing produces worse sensitivities than not summing if the duty cycles 
of the pulsations are large. This results from the fact that such pulsations have only a small number of significant harmonics 
(see Fig. 4), so that summing tends to add only noise rather than signal. 
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Fourier Frequenoy — 566000 Power / Local Power 

Fig. 6. — An 18 <T (single trial) detection from the discovery observation of the 1.7 hour binary PSR J1807— 2459 in the globular 
cluster NGC6544 using a Fourier-domain "acceleration" search. Contour intervals correspond to 30, 60, 90, 120, and 150 times 
the average local power level. The intrinsic pulsar period and / = (which corresponds to an un-accelerated FFT of the data) 
are marked by the solid gray lines. The dots correspond to the "raw" or un-interpolatcd powers from the original FFT of the 
observation. The gray ellipse is the predicted "path" of the pulsar in the /-/ plane given the known binary parameters. During 
the 28.9 minute observation, the pulsar moved from ~ 11 o'clock to ~ 3 o'clock on the ellipse. The peak's slight offset from 
the ellipse and the presence of "shoulders" indicate that the constant / assumption of the acceleration search could not fully 
correct for the orbital motion during this observation. The top and right-hand panels show cuts through the peak in the / 
and / directions respectively. The line in the top panel with gray dots shows the Fourier interpolated / = power spectrum 
(calculated as per §4.1). 
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Fig. 7. — The theoretical response for a boxcar-windowed signal with a constant frequency derivative. Contours plotted are 
10%, 30%, 50%, 70%, and 90% of peak (i.e. fully corrected) power. The top panel shows the familiar sine response of the signal 
along the / = line. The right panel shows a similar cut along the Ar = line (i.e. the calculated average frequency is the 
true average signal frequency) . The relatively uniform spread of signal power over the local power spectrum is apparent for all 
values of frequency derivative. 



